######################################################################################
##Replication Code for ###############################################################
##Decision by Design: Leaders, Bureaucracies, and International Crisis Performance ###
##Tyler Jost #########################################################################
##International Studies Quarterly ####################################################
##(2) Analysis #######################################################################
######################################################################################

rm(list = ls())
require(here)
require(interplot)
require(sandwich)
require(clubSandwich)
require(lmtest)
require(stargazer)
require(xtable)
require(ggplot2)
require(ggpubr)

#########################################################################
##Figure 1 ##############################################################
#########################################################################

rm(list = ls())
load("figure_1.RData")

##Regime Type - Autocracy
holder <- as.data.frame(table(d.descriptive$type, d.descriptive$gwf_regimetype))
autocracy <- holder[holder$Var2!="Democracy",]
p1 <- ggplot(autocracy, aes(x=Var2, y=Freq, fill=Var1)) +
  geom_bar(stat="identity", position=position_dodge()) + ylim(0,850) + scale_fill_brewer(palette = "Blues") +
  theme_minimal() + xlab("") + ylab("Number of Country-Years") + ggtitle("Autocracy") + theme(plot.title = element_text(hjust = 0.5, face="bold"), legend.position = "none", legend.title = element_blank())


##Regime Type - Democracy
democracy <- as.data.frame(table(d.descriptive$type[d.descriptive$regime!="Dictatorship"], d.descriptive$regime[d.descriptive$regime!="Dictatorship"]))
p2 <- ggplot(democracy, aes(x=Var2, y=Freq, fill=Var1)) +
  geom_bar(stat="identity", position=position_dodge()) + ylim(0,850) + scale_fill_brewer(palette = "Blues") +
  theme_minimal() + xlab("") + ylab("") + ggtitle("Democracy") + theme(plot.title = element_text(hjust = 0.5, face="bold"), legend.position = "none", legend.title = element_blank())


combined <- ggarrange(p1, p2, ncol = 2, nrow = 1, common.legend = TRUE, legend = "bottom")

ggsave("figure_1.pdf", plot = combined, scale = 1, width = 10, height = 3.5, dpi = 500, limitsize = T)

#########################################################################
##Load Crisis Data ######################################################
#########################################################################

rm(list = ls())
load("analysis.RData")

#########################################################################
##Descriptive Statistics on Crises ######################################
#########################################################################

#Number of unique crisis initiations
nrow(d[d$crisis_initA==1,]) #683 observations

#Number of unique crises
length(unique(d$crisno)) #343 crises total

#Number of multilateral crises
multiple <- table(d$crisno[d$crisis_initA==1])
length(multiple[multiple>1]) #132 crises involve more than one initiator

#Cross-tabs
cross <- table(d$crisis_init_failA[d$crisis_initA==1], d$type[d$crisis_initA==1])
t1 <- xtable(cross)
t2 <- xtable(round(prop.table(cross, 2), 2))

#########################################################################
##Table 1 ###############################################################
#########################################################################

R2logit <- function(y,model){
  R2<- 1-(model$deviance/model$null.deviance)
  return(R2)
}

crse <- function(x, df){
  d.cse <- na.omit(df[, c("crisno", all.vars(formula(x)))])
  fc <- d.cse$crisno
  m <- length(unique(fc))
  k <- length(coef(x))
  u <- estfun(x)
  u.clust <- matrix(NA, nrow=m, ncol=k) 
  for(j in 1:k){
    u.clust[,j] <- tapply(u[,j], fc, sum)
  }
  bread <-vcov(x)
  c.vcov <- bread %*% ((m / (m-1)) * t(u.clust) %*% (u.clust)) %*% bread
  result <<- coeftest(x, c.vcov)
  return(result)
}

##Model 1: 
m1.1 <- glm(crisis_init_failA ~ type + as.factor(statea), 
                   family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR1.1 <- round(R2logit(crisis_init_failA, m1.1),3)
obs1.1 <- nobs(m1.1)
m1.1 <- crse(m1.1, d[d$crisis_initA==1,])

##Model 2:
m1.2 <- glm(crisis_init_failA ~ type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
                    family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR1.2 <- round(R2logit(crisis_init_failA, m1.2),3)
obs1.2 <- nobs(m1.2)
m1.2 <- crse(m1.2, d[d$crisis_initA==1,])

##Model 3:
m1.3 <- glm(crisis_init_failA ~ type + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + size + I(size^2) + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR1.3 <- round(R2logit(crisis_init_failA, m1.3),3)
obs1.3 <- nobs(m1.3)
m1.3 <- crse(m1.3, d[d$crisis_initA==1,])

stargazer(m1.1, m1.2, m1.3, type = "latex", omit=c("Constant", "time", "time2", "time3", "statea"),
          omit.stat=c("f", "ser", "ll", "aic"),
          label=c("table: main_model"),
          covariate.labels = c("Siloed", "Fragmented", "Dictatorial", "Leader Age", "Male", "Military Experience", "Combat Experience", "Time in Office", "Military Regime",
                               "Personalist Regime", "Regime Transition", "State A Democraticness", "State B Democraticness", "State A Capabilities", "State B Capabilities",
                               "Relative Capabilities", "Alliance Score", "Body Size", "Body Size -- Quadratic"),
          add.lines=list(c("Observations", obs1.1, obs1.2, obs1.3),
                         c("McFadden Pseudo-R$^2$", pseudoR1.1, pseudoR1.2, pseudoR1.3),
                         c("Cluster Robust SEs", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Country Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$")),
          dep.var.labels = c("Failure to Achieve Goals"),
          title = "National Security Institutions and International Crisis Performance",
          notes.append=FALSE, no.space=T, single.row = T,
          notes = c("Robust standard errors are clustered by crisis.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))

#########################################################################
##Check differences between institutional types #########################
#########################################################################

##Number of dictatorial crisis intiations
round(nrow(d[d$type=="Dictatorial" & d$crisis_initA==1,]) / nrow(d[d$crisis_initA==1,]),3) #6\%

##Check difference between siloed and fragmented
d$type2 <- relevel(d$type, ref="Siloed") 

##Siloed as base: 
m4 <- glm(crisis_init_failA ~ type2 + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
crse(m4, d[d$crisis_initA==1,]) # p < 0.002

#########################################################################
##Continuous Measure ####################################################
#########################################################################

##Model 1: 
m2.1 <- glm(crisis_init_failA ~ info_up + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR2.1 <- round(R2logit(crisis_init_failA, m2.1),3)
obs2.1 <- nobs(m2.1)
m2.1 <- crse(m2.1, d[d$crisis_initA==1,])

##Model 2: 
m2.2 <- glm(crisis_init_failA ~ info_horizontal + as.factor(statea), 
            family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR2.2 <- round(R2logit(crisis_init_failA, m2.2),3)
obs2.2 <- nobs(m2.2)
m2.2 <- crse(m2.2, d[d$crisis_initA==1,])

##Model 3: 
m2.3 <- glm(crisis_init_failA ~ info_up + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
          family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR2.3 <- round(R2logit(crisis_init_failA, m2.3),3)
obs2.3 <- nobs(m2.3)
m2.3 <- crse(m2.3, d[d$crisis_initA==1,])

##Model 4: 
m2.4 <- glm(crisis_init_failA ~ info_horizontal + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
            family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR2.4 <- round(R2logit(crisis_init_failA, m2.4),3)
obs2.4 <- nobs(m2.4)
m2.4 <- crse(m2.4, d[d$crisis_initA==1,])

##Model 5: 
m2.5 <- glm(crisis_init_failA ~ info_up + info_horizontal + age + gender + milservice + combat + hos_experience + gwf_military + gwf_personal + regtransA + polity2A + polity2B + LcincA + LcincB + cinc_rel + s_wt_atop + as.factor(statea), 
            family=binomial(link = "logit"), data = d[d$crisis_initA==1,])
pseudoR2.5 <- round(R2logit(crisis_init_failA, m2.5),3)
obs2.5 <- nobs(m2.5)
m2.5 <- crse(m2.5, d[d$crisis_initA==1,])

stargazer(m2.1, m2.2, m2.3, m2.4, m2.5, type = "latex", omit=c("Constant", "statea"),
          omit.stat=c("f", "ser", "ll", "aic"),
          label=c("table: continuous_measures"),
          covariate.labels = c("Information Search", "Bureaucratic Access", "Leader Age", "Male", "Military Experience", "Combat Experience", "Time in Office", "Military Regime",
                               "Personalist Regime", "Regime Transition", "State A Democraticness", "State B Democraticness", "State A Capabilities", "State B Capabilities",
                               "Relative Capabilities", "Alliance Score"),
          add.lines=list(c("Observations", obs2.1, obs2.2, obs2.3, obs2.4, obs2.5),
                         c("McFadden Pseudo-R$^2$", pseudoR2.1, pseudoR2.2, pseudoR2.3, pseudoR2.4, pseudoR2.5),
                         c("Cluster Robust SEs", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"),
                         c("Country Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$")),
          dep.var.labels = c("Failure to Achieve Goals"),
          title = "National Security Institutions and International Crisis Performance (Continuous Measure)",
          notes.append=FALSE, no.space=T, single.row = T,
          notes = c("Robust standard errors are clustered by crisis.$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}"))

